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Abstract 

We present a formulation of the Constrained Path Monte Carlo (CPMC) 
method for fermions that uses trial wave-functions that include many-body 
effects. This new formulation allows us to implement a whole family of gen- 
eralized mean-field states as constraints. As an example, we calculated su- 
perconducting pairing correlation functions for the two-dimensional repulsive 
Hubbard model using a BCS trial state as the constraint. We compared the 
results with the case where a free-electron trial wave- function is used. We 
found that the correlation functions are independent of which state is used as 
the constraint, which reaffirms the results previously found by Zhang et. aJU 
regarding the suppression of long range pairing correlations as the system size 
increases. 

PACS Numbers: 74.20.-z, 74.10.+V, 71.10.Fd, 71.10.-w, 02.70.Lq 
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I. INTRODUCTION 



Since the discovery of high temperature superconductivity, an enormous effort has been 
devoted to the theoretical study of two-dimensional electronic models. This effort is driven 
by the belief that the mechanism for superconductivity lies within the CuC>2 planes common 
to these materials and is dominantly electronic in origin. The two-dimensional repulsive 
Hubbard model has attracted the most attention as the simplest effective model possibly 
embodying the key electronic phenomena at low energies. Numerous works on this model 
have reproduced qualitatively the observed magnetic properties of the cuprates in the normal 
state.! However, the search for superconductivity in the Hubbard model, although intensive 
and extensive, has yielded few positive indicators.! 

Most of the present knowledge on the phase diagram of the two-dimensional repulsive 
Hubbard model has been obtained by combination of theorems and numerical studies of 
finite size clusters. The numerical studies used Lanczos, Variational Monte Carlo, and zero 
or finite temperature quantum Monte Carlo techniques. In a superconducting phase, one 
expects the superconducting pairing correlation functions to exhibit off-diagonal long range 
order (ODLRO), which is an indication of the Meissner effect.! With this in mind, a number 
of investigators have calculated pairing correlation functions in various symmetry channels. 
However, most calculations were limited to high temperatures and small system sizes. In 
the case of Monte Carlo studies these limitations were imposed by the fermion sign problem 
which causes the variances of computed quantities and hence the computing time to grow 
exponentially with the increase in system sizes. 

Recently, a new zero temperature quantum Monte Carlo method, the Constrained Path 
Monte Carlo (CPMC), was developed that overcomes the major limitations of the sign 
problem.! This method allows the calculation of pairing correlation functions at zero tem- 
perature without the exponential increase in computer time with system size. Using this 
method, Zhang et al.l calculated <ij.2_.y2 -wave and extended s-wave pairing correlation func- 
tions versus distance in the ground state for lattices up to 16 x 16. They found that the 
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d x 2_ y 2-wdNe correlations are stronger than extended s-wave correlations. However, as the 
system size or the interaction strength was increased, the magnitude of the long-range part 
of both correlation functions vanished. 

Although the findings of Zhang et al.S provide evidence for the absence of ODLRO in the 
two-dimensional Hubbard model, the CPMC method is approximate and has a systematic 
error which is difficult to gauge. The systematic error is associated with the wave-function 
used to constrain the Markov chains produced by the Monte Carlo procedure. More specif- 
ically, in the CPMC method the ground state wave-function is represented by an ensemble 
of Slater determinants. As these determinants evolve in imaginary time, the ones with a 
negative overlap with a constraining wave-function are discarded. This procedure elimi- 
nates the sign problem but introduces an approximation that depends on the quality of the 
constraining wave-function. Zhang et al.0 used free-electron and unrestricted Hartree-Fock 
wave-functions. More sophisticated choices of wave-functions, particularly ones exhibiting 
strongly correlated electron effects, are typically difficult to implement, because of the in- 
creasing number of Slater determinants needed and the consequent increase in computing 
time. 

In this work, we extended the formulation of the CPMC method in a way that allows 
the use of a wide variety of trial wave-functions with only a small increase in computing 
time. As an illustration, we calculated the superconducting pairing correlation functions of 
the two-dimensional repulsive Hubbard model in the rf x 2_ 2/ 2-wave channel using as a con- 
strain a BCS wave-function that has superconducting ODLRO. We found that the resulting 
correlation functions are the same as those obtained using the free-electron and Hartree- 
Fock constraining wave-functions. This reaffirms the results by Zhang et al.0 regarding the 
vanishing of long range pairing correlations as the system size increases. 

The article is organized as follows: in section [Tl] we briefly describe the CPMC technique 
emphasizing aspects of the new formulation. In section [TTT] we define the Hamiltonian and 
pairing correlation functions and present our results. In section [TV] we discuss our conclu- 
sions. 
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II. METHOD 



In this section we summarize the main features of the CPMC method. For a more detailed 
description of the method see Ref. |j. In the CPMC method, the ground-state wave-function 
l^o) is projected in imaginary time r from a known initial wave-function \^f(r — 0)) — \^t) 
by a branching random walk in an over-complete space of Slater determinants \<p), 

N a N 

!</>> = 11410) ; 0l = £4<&Ji> (l) 

i,a j=l 

where cj CT creates and electron in orbital j with spin a (rij a = Cj^Cj^), and 

m 1 ) * <w (2) 

with iV the number of available single-particle states (for the Hubbard model corresponds 
to the total number of lattice sites) and N a the number of particles with spin a. The total 
number of electrons is given by iV e = iVf + iVj_. 

The projection corresponds to finding the ground-state from the long-time solution of 
the imaginary-time representation of Schrodinger's equation specified by a Hamiltonian H 

^± = -(H-E l)\V) (3) 

with E the ground-state energy (h is set to 1). 

Provided No = (^ol^O)) 7^ and H being time-independent, the formal solution 

|tf( r )) = e- T ^~ Eoi) \y(0)} (4) 

has the property 

lim |*(r)> = iV |^ ) (5) 

T — >00 

On the computer this large r limit is accomplished by breaking up r in small time-steps At 
and iterating the equation 



where Et is a guess at the ground-state energy E$ and AtN s = r with iV s the number of 
imaginary time-steps. As r —>■ oo, the iteration becomes stationary, i.e. d\^)/dr = 0, and 
if Et is adjusted to equal E , then |^(r — > oo)) = A^l^o)- 

The propagation in imaginary time is done in the following way: in the space of Slater 
determinants, we write \^ Q ) = X^x(0)|0) an d choose > 0. By being positive, the 

function describes the distribution of Slater determinants representing the ground state. 
The Monte Carlo process samples from this distribution. This process is implemented by 
the application of a Trotter decomposition and a Hubbard-Stratonovich transformation to 
the iterative equation (^j) and converting it into 

|^ n+1 ) = J dxP(x) J B(x)|^ n ) (7) 

where x is a multi-dimensional random variable distributed according to -P(x) and -B(x) is 
an operator approximating e~ ArH for a given value of the random variable, whose general 
structure is a product of exponentials of operators quadratic in c and c'. For each time step 
At, -B(x) has the property of transforming one Slater determinant into another. The Monte 
Carlo method evaluates the mult i- dimensional integral (|7|) by using an ensemble of random 
walkers represented by Slater determinants \<p). For each walker, it samples x from -P(x) 
and then obtains the new Slater determinant by multiplying 

\4> n+l ) = B{*W) (8) 

Once the Monte Carlo procedure converges, the ensemble of |0) represents \^o) in the sense 
that their distribution is a Monte Carlo sampling of in this sense, the CPMC approach 
is a sort of stochastic configuration interaction method. 

To specify the ground-state wave-function completely, only determinants satisfying 
(\&o|0) > are needed because |\l/o) resides in either of two degenerate halves of the Slater 
determinantal space (in general, a manifold of dimension N e (N — N e )) , separated by a nodal 
hypersurface M defined by (^o|0) = 0. The sign problem occurs because walkers can cross 
M as their orbitals evolve continuously in the random walk. Asymptotically in r they pop- 
ulate the two halves equally, leading to an ensemble that tends to have zero overlap with 
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I^o)- If were known, one would simply constrain the random walk to one half of the space 
and obtain an exact solution of Schrddinger's equation. In the CPMC method, without a 
priori knowledge of A/", we use a constraining wave-function, which we usually take to be 
the trial wave-function \^t), and require the Slater determinants to satisfy > 0. 

Thus, the quality of the calculation clearly depends on \^t)- in the past only free-electron 
or Hartree-Fock wave-functions were implemented, mainly due to their simplicity and the 
novelty of the method. However, it is desirable to use more sophisticated wave-functions 
that include many-body effects. For example, to study superconductivity it is interesting to 
implement trial wave-functions that exhibit ODLRO, like a BCS wave-function. 

Our goal is to use trial wave-functions of the type (i. e., a Bogoliugov transformation of 
the vacuum |0), (0|0) = 1) 

I*t> = IIK + ^4A|)|0> (9) 
k 

where the product includes all values of momentum k = (k x , k y ) in the first Brillouin zone 
and \uk\ 2 + \vk\ 2 = 1 to ensure normalization ((^t\^t) — !)• Other than satisfying the 
normalization condition, the parameters u k and Vk can be chosen arbitrarily. 

Equation (^) represents a wave-function that does not have a fixed particle number 
N e . To represent a fixed electron number, \^t) needs to be projected onto that particular 
subspace. The resulting wave-function is a linear combination of a large number of Slater 
determinantsi (large in the sense that the number grows very rapidly with system size and 
particle number to the point where it becomes impractical to use). Alternatively, one can 
work in an extended space with different electron numbers. To do that, we follow Yokoyama 
and ShibJ and perform a particle-hole transformation on one of the spin species: 

dk = cl kl 

ft 

C k — C k\ 

Using this transformation and noting that the new vacuum |0) is related to the old one by 

io)=n4io) (ii) 



we can rewrite \^t) in terms of the new c and d operators: 

i*T>=n(«*4+«*4)io> (12) 

k 

so that \^t) is represented by a single Slater determinant. Since we are interested in 
projecting out the ground state with a fixed electron number, we have to use the propagator 
e - T ( H -Eot-nN e ) _ £Y( T ) anc i choose /x, the chemical potential, to select the desired number 
of electrons N e = (^ |iV e |^ )/(^ |^ ) (N e = J2ja n ja)- At the end of the projection the 
ground state wave-function will have a fixed number of electrons given by the choice of fi. 

The changes in the CPMC method necessary to use the BCS form of a correlated wave- 
function are minor. Instead of matrices $°" for up and down spin of sizes N x N a to represent 
the random walkers, they, as well as the trial wave- function \^t), are now represented by a 
single matrix of size 2N x N. The increase in computation time caused by the increase in 
the size of the matrices depends on the system size and the number of particles. A rough 
estimate gives the increase as the factor 3N/N e . For example, for a 6 x 6 system with 
N e = 26 this is 4 = 2.89N/N e . The closer we get to half-filling (JV e = N) the smaller the 
increase. In general, for the filling fractions studied here, the increase in computer time is 
of the order of 4. 



III. CALCULATION AND RESULTS 

The Hamiltonian is the usual Hubbard Hamiltonian on a square lattice with periodic 
boundary conditions: 

# = -* £ (cl^ + c^c^ + U^nqn^ (13) 

<ij>,tr i 

where t is the nearest neighbor hopping matrix element and U is the on-site Coulomb 
repulsion. We set t — 1 so that all energies are measured in units of t. In terms of the 
operators c and d defined by the transformation (|i~0|) the Hamiltonian has the form 

H = -t £ (cjc. + cjq - d\d 3 - did,) + U £ <(1 - nf) (14) 

<ij> i 
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where n\ (nf) denotes the occupation in the c (d) orbital. This transformed Hamiltonian 
corresponds to a two-band spinless fermion model. 

We computed the ground-state energy and the superconducting pairing correlation func- 
tions in the d x 2_ y 2 -wave channel using the following definitions: 

P d (R) = (Al(R)A d (0)) (15) 

where the pair field operator is 

A d (R) = E fd{S) [c^c^ - c^c^] (16) 

6 

with 5 = ±x, ±y , fd(±x) = 1 and fd{±y) — —1 . R denotes the position in the lattice in 
units of the lattice constant which is taken to be unity. 

We used trial wave-functions of the form (|9|) with u k and v k given by the BCS relation 

Vk = A fc 

u * e k -fi+^/(e k ~fi) 2 + \A k \ 2 

where e k is a single particle energy and A^ is the gap, = Af(k). A is a variational 
c-number and f(k) represents the symmetry of the pairing which we choose to be d x 2_ y 2, 
f(k) = cos(k x ) - cos(ky). 

We concentrated in the d^.^-wave channel in part because the existence of ODLRO in 
the extended s-wave channel is conditioned upon the existence of ODLRO in the isotropic 
s-wave channel.S Since the possibility of pairing in the isotropic s-wave channel is highly 
unlikely for the repulsive Hubbard model, so is the chance of pairing in the extended s-wave 
channel. Moreover, these statements have been verified numerically by us and by Zhang et 
al.0 Also, it has been increasingly established experimentally that the order parameter in 
the superconducting cuprates has fi^.^-wave symmetry. 

We used two different trial wave-functions: one with A = 0.5, which corresponds to 
a BCS superconducting state, and the other one with A = 0, which corresponds to the 
free-electron case. In both cases we choose the parameter \l in the BCS wave-function so 
that (^rl-^Vel^r) = -^e where N e is the number of electrons we are interested in. While the 
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free-electron wave-function has a fixed number of electrons (cr/v e = \J (N e ) — (N e ) 2 = 0), 
the BCS wave-function with A 7^ has components with different electron numbers so that 
a N e 7^ 0. It is important to notice that in general the parameter fi in the BCS wave-function 
is different than the one used in the propagator hi{r). The latter one is set so that at the 
end of the propagation the ground state has the desired number of electrons N e . 

To illustrate the difference between these two wave-functions, in Fig. 1 we plot 
the variational value of the d^^-wave correlation functions versus distance, that is 
( 1 I r y|A^(i2)A £ i(0)| 1 J r r), for the two trial wave-functions in a 10 x 10 system with U = 4 
and N e = 82, so that the filling fraction is n e = N e /N = 0.82. This filling corresponds to a 
closed shell case, that is, the free-electron ground state is non-degenerate. In the free-electron 
case the correlations die out rapidly with distance, while in the BCS case the existence of 
ODLRO is evident in the sense that for long distances, the correlation functions approach 
a finite value given by the square of the superconducting order parameter A : 

a 5C = I E f(k> k v k = ± £ /(*) T . \ A2 (is) 

The overlap between the two normalized trial wave-functions is (^r(A = 0)|^t(A = 0.5)) = 
0.0076, so the two wave-functions are close to being orthogonal. 

The variational energy E v = (^/t\H\^t) is much larger for the BCS trial wave-function 
than for the free-electron trial wave-function. In general we find that the variational energy 
increases monotonically with the parameter A of the BCS wave-function, as it is shown in 
Fig. 2 for a 10 x 10 system with U = 4 and (N e ) = 82. This variation contrasts previous 
results obtained with the Variational Monte Carlo method, which found that a non-zero 
value of A minimizes the variational energy.Stl'EI However, in these Gutzwiller factor 

was included in the wave-function that projected out totally or partially the states with 
double occupancy. It seems that the inclusion of this factor is crucial to obtain a minimum 
of the variational energy at a finite value of A. At present, our formulation does not allow 
the use of trial wave-functions that are non-Fock states such as the Guztwiller wave-function: 
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\^ G ) = \[(i-gn^n i[ )\m F ocK) (19) 

i 

with g a variational parameter that determines the average number of doubly occupied 
sites. (When g = 1, double occupation is completely suppressed.) Even though such wave- 
functions are not implemented, since we are doing a projection in imaginary time onto the 
ground state of the system, it is not crucial to improve the variational energy of our trial 
state. 

In the large U limit, the Hubbard model can be mapped onto the t — J model . This 
strong coupling limit was used in Refs. |8] and |6| to calculate the energy, making a comparison 
with our work difficult. However, we can do a comparison with Ref. ^| since they used the 
Hubbard Hamiltonian to calculate the energy. In their Fig. 1 they report the variational 
energy per site as a function of A for a 6 x 6 system with U = 8, 32 electrons, periodic 
boundary conditions in the x direction and anti-periodic in the y direction. From their 
figure, the minimum value for the energy per site is -0.65523 and corresponds to a value 
of A = 0.1. The variational energy per site that we obtain for the same system but with 
periodic boundary conditions in both directions is 0.02726. The difference can likely be 
accounted for by the fact that we did not project our wave-function onto a fixed particle 
number and second, we did not use a Gutzwiller factor. However, the ground state energy 
per site calculated with the CPMC method is —0.7272 ±0.0005, which is considerably lower 
than their value. 

As a check of our algorithm we compared the correlation functions and ground-state 
energy given by the CPMC method using the free-electron trial wave-function with results 
by Zhang et al.,0 who used the original formulation of the CPMC, for a 6 x 6 system with 
U = 4 and N e = 26 and an 8 x 8 system with U = 8 and N e = 50. We found excellent 
agreement with their results. 

In Fig. 3 we plot the resulting correlations functions given by the CPMC calculation 
with the two trial wave-functions used in Fig. 1, for 10 x 10 with U — 4. It is clear that 
the results are essentially the same no matter what trial wave-function is used. The long 



10 



distance magnitude of the correlation functions is very small, smaller than the free-electron 
case. 

Similar calculations to the ones presented in Fig. 3 were done for 8 x 8 and 6x6 
systems with U — 4, 6 and 8 and dopings corresponding to closed shells cases. The results 
are consistently the same: the correlation functions are the same no matter what trial 
wave-function is used. The ground-state energy, however, is always larger when the BCS 
wave-function is used. The difference between the two ground-state energies is larger for 
larger U. When the BCS wave-function is used, we find that there are more nodal crossings; 
that is, more walkers are discarded because their overlap with the trial wave-function is 
negative. We believe this is why the energy is higher in the case of the BCS wave-function. 

We did not use systems larger than 10 x 10 in part because as system size increases, it 
becomes more difficult to select ji in the propagator to get the desired number of electrons. 
This is because the energy levels are getting closer in larger systems. Also, we found that 
the correlation functions are the same no matter which trial wave-function is used for 6 x 6, 
8x8 and 10 x 10 systems. This evidence is enough to conclude that the correlation functions 
are independent of which trial wave-functions is used. 

IV. CONCLUSIONS 

We presented a formulation of the CPMC method that uses trial wave-functions that 
include correlation effects and have components of different electron numbers. Instead of 
projecting it onto a subspace with fixed number of electrons, we used a particle-hole trans- 
formation in one of the spin species to write such trial wave-functions as only one Slater 
determinant. 

Because of the increase in the size of the matrices used, this formulation involves a small 
increase in computing time compared to the original formulation. The increase in CPU 
time is roughly 3N/N e . For the dopings considered in this work it comes to a factor of 
approximately 4. 
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This new formulation is very general and allows the implementation of a whole family of 
mean-field wave-functions. Following Bach, Lieb and Solovej0 we call this class of functions 
generalized Hartree-Fock states, i. e., states that are ground states of some quadratic mean- 
field Hamiltonian in Fock space which do not necessarily conserve particle number. Possible 
examples include spin-density wave, charge-density wave and superconductivity. 

As an illustration, and because of its importance in high temperature superconductivity, 
we used a BCS trial wave-function with d x 2 _ y 2 -wave symmetry to calculate the supercon- 
ducting pairing correlation functions in the ground state for the two-dimensional repulsive 
Hubbard model. We compared this result with the one using the free-electron trial wave- 
function. We studied 6 x 6, 8 x 8, and 10 x 10 systems for different values of U and dopings 
and found that the results for the correlation functions are independent of which trial wave- 
function is used for the constraint. 

Most of the calculations presented in this work correspond to closed shell cases, that is, 
electron fillings with a non-degenerate free-electron ground state. To check the consistency 
of our results we also studied some open shell cases like a 6 x 6 system with 32 electrons 
(n e = 0.89), U = 8 and periodic boundary conditions. We used three different trial wave 
functions: one free-electron wave function with a fixed number of electrons, another free- 
electron wave function but with some paired electrons in the Fermi surface and a BCS 
wave-function with A = 0.1. The CPMC result is consistent with those of the closed shell 
cases: the superconducting pairing correlation functions, which vanish for large distances, 
are independent of the trial wave-function used. Technically, the open shell case is more 
difficult because in general the free-electron trial wave-functions do not have translational 
invariance. For this reason, one finds different values of the correlation functions for the same 
distance \R\ but different directions in the lattice. To overcome this problem we averaged 
the correlation functions for a given \R\ over all possible directions in the lattice. This 
procedure is also used for the closed shell cases but is more relevant in the open shell case 
where the differences are caused by a broken symmetry introduced by the trial wave-function 
as oppossed to small statistical fluctuations due to the Monte Carlo process. 

12 



These results reaffirm the previous ones by Zhang et al.EJ implying the absence of ODLRO 
in the cL^^-wave channel of the two-dimensional repulsive Hubbard model. We do not 
dismiss the possibility of ODLRO existing in some exotic channel or for some combination 
of quasiparticle operators instead of the bare ones.0 This work has only investigated the 
channels commonly studied. Although it is not rigorously proven that the absence of ODLRO 
implies no Meissner effect and consequently no superconductivity, it is reasonable to think 
that a model without apparent ODLRO is inappropriate as a model of the superconducting 
phase for the high temperature superconducting materials. 

The lack of clear numerical evidence of d x 2_ y 2-wave superconductivity upon doping and 
the abundance of clear numerical evidence of antiferromagnetism at half filling makes it 
hard to see how a theory, like the SO (5) phenomenology, can apply to the Hubbard model 
as some have recently suggested.^ This phenomenology requires the antiferromagnetic long 
range order at half-filling to transform into d^.^-wave superconducting long range order 
in the doped states. If the low lying excited states have approximate SO (5) symmetry, why 
then does the strong antiferromagnetic state transform into something that is so hard to 
find? The two-dimensional repulsive Hubbard model seems to be an inappropriate candidate 
for the SO (5) phenomenology. 

V. ACKNOWLEDGMENTS 

We are thankfull to S. Trugman for the critical reading of the manuscript. The C++ 
program used for this work incorporated the MatrixRef matrix classes written by S.R. White, 
available at [http:/ /hedrock.ps. uci.edu} This work was supported by the Department of 
Energy. Some of the calculations were performed on the computers at NERSC. 



13 



REFERENCES 

1 Shiwei Zhang, J. Carlson and J. E. Gubernatis, Phys. Rev. Lett. 78, 4486 (1997). 

2 See for example, Elbio Dagotto, Rev. Mod. Phys. 66, 763 (1994). 
3 G. L. Sewell, J. Math. Phys. 38, 2053 (1997). 

4 Shiwei Zhang, J. Carlson and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995); Phys 
Rev. B 55, 7464 (1997). 

5 J. P. Bouchaud, A. Georges and C. Lhuillier, J. Phys. France 49, 553 (1988). 
6 H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 57, 2482 (1988). 
7 Shoucheng Zhang, Phys. Rev. B 42, 1012 (1990). 

8 Claudio Gross, Phys. Rev. B 38, 931 (1988). 

9 T. Nakanishi, K. Yamaji, T. Yanagisawa, J. Phys. Soc. Jpn. 66, 294 (1997). 

10 Shiwei Zhang, J. Carlson and J. E. Gubernatis, private communication. 

11 V. Bach, E. H. Lieb and J. P. Solovej, J. Stat. Phys. 76, 3 (1994). 

12 W. Hanke, R. Eder, E. Arrigoni, A. Dorneich, S. Meixner and M. G. Zacher, |cond- 
mat/98070T5 . 



13 



E. Dagotto and J. R. Schrieffer, Phys. Rev. B 43, 8705 (1991). 



14 



FIGURES 

FIG. 1. Variational value of the pairing correlations versus distance \R\ for two different trial 
wave-functions in a 10 x 10 system. Parameters are U = 4 and filling fraction n e = 0.82. The BCS 
wave- function exhibits ODLRO. 

FIG. 2. BCS variational energy per site as a function of A for the same system as in Fig. 1. 
The energy increases monotonically with A. The inset shows smaller values of A where Ref. 9 
finds a minimum. 

FIG. 3. Pairing correlation functions in the d a ,2_ y 2-wave channel given by the CPMC method 
for same system as in Fig. 1. The inset shows the long range part in detail. The results are the 
same for the two different trial wave-functions: the correlations decay quickly with distance. Errors 
bars are smaller than the size of the symbols. 
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